%% 原始信号
% 定义.mp3文件的路径
mp3FilePath = './music.mp3';
% 调用getMusic函数加载音频数据
[audioData, sampleRate] = getMusic(mp3FilePath); % 确保这里是mp3FilePath
% 打印采样频率和采样点数
fprintf('采样频率: %d Hz\n', sampleRate);
fprintf('采样点数: %d\n', length(audioData));
% 绘制原始音频波形图
t = (0:length(audioData)-1) / sampleRate; % 创建时间向量
figure; % 创建一个新的图形窗口
plot(t, audioData);
xlabel('时间 (秒)');
ylabel('幅度');
title('原始音频波形图');
% 进行FFT变换
Y = fft(audioData);
% 计算双边频谱 P2
P2 = abs(Y/length(audioData));
% 计算单边频谱 P1
halfIndex = floor(length(audioData)/2); % 使用floor函数确保索引是整数
P1 = P2(1:halfIndex+1);
P1(2:end-1) = 2*P1(2:end-1);
% 计算频率轴
f = sampleRate*(0:halfIndex)/length(audioData);
% 绘制原始音频频谱图
figure;
plot(f, P1);
xlabel('频率 (Hz)');
ylabel('幅度');
title('原始音频频谱图');
% 定义滤波器的参数
Rp = 1;    % 通带波动（dB）
Rs = 100;  % 阻带衰减（dB）

%% FIR低通
% 低通滤波器
fp = 1000; % 通带截止频率（Hz）
fs = 1200; % 阻带截止频率（Hz）
% 归一化频率
Wp = 2*fp/sampleRate;
Ws = 2*fs/sampleRate;
% 计算过渡带宽
DeltaW = Ws - Wp;
% 使用Hamming窗的经验公式计算滤波器阶数N
N_low = ceil((6.6*pi) / DeltaW);
% 设计FIR低通滤波器
b_low = fir1(N_low, Wp, 'low', hamming(N_low+1));

%% 评估
% 计算频率响应
[H, f0] = freqz(b_low, 1, 1024, sampleRate);
% 绘制幅度响应
figure;
subplot(2, 1, 1);
plot(f0, abs(H));
title('FIR低通滤波器的幅度响应');
xlabel('频率 (Hz)');
ylabel('幅度');
% 绘制相位响应
subplot(2, 1, 2);
plot(f0, angle(H));
title('FIR低通滤波器的相位响应');
xlabel('频率 (Hz)');
ylabel('相位 (弧度)');
% 绘制群时延
groupDelay = grpdelay(b_low, 1, 1024, sampleRate);
figure;
plot(f0, groupDelay);
title('FIR低通滤波器的群时延');
xlabel('频率 (Hz)');
ylabel('群时延 (样本)');
% 绘制滤波器的冲激响应
impulseResponse = impz(b_low, 1);
figure;
stem(impulseResponse);
title('FIR低通滤波器的冲激响应');
xlabel('样本');
ylabel('幅度');

% 对音频信号进行低通滤波处理
filteredAudioData_low = filter(b_low, 1, audioData);
% 绘制低通滤波后的音频波形图
figure;
plot(t, filteredAudioData_low);
xlabel('时间 (秒)');
ylabel('幅度');
title('FIR低通滤波后的音频波形图');
% 对低通滤波后的音频信号进行FFT变换
Y_low = fft(filteredAudioData_low);
% 计算双边频谱 P2_low
P2_low = abs(Y_low/length(filteredAudioData_low));
% 计算单边频谱 P1_low
P1_low = P2_low(1:length(P2_low)/2+1);
P1_low(2:end-1) = 2*P1_low(2:end-1);
% 绘制低通滤波后的音频频谱图
figure;
plot(f, P1_low);
xlabel('频率 (Hz)');
ylabel('幅度');
title('FIR低通滤波后的音频频谱图');
%% butter低通
% 计算滤波器阶数
[N, Wn] = buttord(Wp, Ws, Rp, Rs);
N = min(N, 8); % 限制滤波器阶数
% 使用SOS表示设计滤波器
[sos, g] = butter(N, Wn, 'low');

%% 评估
% 计算频率响应
[H, f0] = freqz(sos, g, 1024, sampleRate);
% 绘制幅度响应
figure;
subplot(2, 1, 1);
plot(f0, abs(H));
title('巴特沃斯低通滤波器的幅度响应');
xlabel('频率 (Hz)');
ylabel('幅度');
% 绘制相位响应
subplot(2, 1, 2);
plot(f0, angle(H));
title('巴特沃斯低通滤波器的相位响应');
xlabel('频率 (Hz)');
ylabel('相位 (弧度)');
% 绘制群时延
groupDelay = grpdelay(sos, g, 1024, sampleRate);
figure;
plot(f0, groupDelay);
title('巴特沃斯低通滤波器的群时延');
xlabel('频率 (Hz)');
ylabel('群时延 (样本)');
% 绘制滤波器的冲激响应
impulseResponse = impz(sos, g);
figure;
stem(impulseResponse);
title('巴特沃斯低通滤波器的冲激响应');
xlabel('样本');
ylabel('幅度');

% 使用filtfilt进行滤波
filteredAudioData_low = filtfilt(sos, g, audioData);
% 绘制滤波后的波形
t = (0:length(audioData)-1)/fs; % 时间向量
figure;
plot(t, filteredAudioData_low);
xlabel('时间 (秒)');
ylabel('幅度');
title('Butterworth低通滤波后的音频波形图');
% 频域分析
Y_low = fft(filteredAudioData_low);
P2_low = abs(Y_low/length(filteredAudioData_low));
P1_low = P2_low(1:length(P2_low)/2+1);
P1_low(2:end-1) = 2*P1_low(2:end-1);
% 绘制频谱图
figure;
plot(f, P1_low);
xlabel('频率 (Hz)');
ylabel('幅度');
title('Butterworth低通滤波后的音频频谱图');

% sound(filteredAudioData_low, sampleRate);  % 播放滤波后的音频
% %% cheby1低通
% % 计算滤波器阶数N
% [N, Wn] = cheb1ord(Wp, Ws, Rp, Rs);
% N = min(N, 4); % 限制滤波器阶数
% % 设计Chebyshev I低通滤波器
% [sos, g] = cheby1(N, Rp, Wn, 'low');
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('cheby1低通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('cheby1低通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('cheby1低通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('cheby1低通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_low = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% t = (0:length(audioData)-1) / sampleRate; % 时间向量
% figure;
% plot(t, filteredAudioData_low);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Chebyshev Type I低通滤波后的音频波形图');
% % 频域分析
% Y_low = fft(filteredAudioData_low);
% P2_low = abs(Y_low/length(filteredAudioData_low));
% P1_low = P2_low(1:length(P2_low)/2+1);
% P1_low(2:end-1) = 2*P1_low(2:end-1);
% % 绘制频谱图
% f = (0:length(P1_low)-1) * sampleRate / length(P1_low); % 频率向量
% figure;
% plot(f, P1_low);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Chebyshev Type I低通滤波后的音频频谱图');
% %% ellip低通
% % 使用ellipord计算最小阶数
% [N, Wn] = ellipord(Wp, Ws, Rp, Rs);
% N = min(N, 4); % 限制滤波器阶数
% % 使用ellip函数设计滤波器
% [sos, g] = ellip(N, Rp, Rs, Wn, 'low'); % 使用SOS表示设计滤波器
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('ellip低通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('ellip低通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('ellip低通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('ellip低通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_low = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% t = (0:length(audioData)-1) / sampleRate; % 时间向量
% figure;
% plot(t, filteredAudioData_low);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Elliptic低通滤波后的音频波形图');
% % 频域分析
% Y_low = fft(filteredAudioData_low);
% P2_low = abs(Y_low / length(filteredAudioData_low));
% P1_low = P2_low(1:length(P2_low) / 2 + 1);
% P1_low(2:end-1) = 2 * P1_low(2:end-1);
% % 绘制频谱图
% f = (0:length(P1_low)-1) * sampleRate / length(P1_low); % 频率向量
% figure;
% plot(f, P1_low);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Elliptic低通滤波后的音频频谱图');
% 
% %% FIR高通
% fs = 4800; % 阻带截止频率（Hz）
% fp = 5000; % 通带截止频率（Hz）
% % 归一化频率
% Wp = 2*fp/sampleRate;
% Ws = 2*fs/sampleRate;
% % 计算过渡带宽
% DeltaW = Wp - Ws;
% % 使用Hamming窗的经验公式计算滤波器阶数N
% N_high = ceil((6.6*pi) / DeltaW);
% % 设计FIR高通滤波器
% b_high = fir1(N_high, Wp, 'high', hamming(N_high+1));
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(b_low, 1, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('FIR高通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('FIR高通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(b_low, 1, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('FIR高通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(b_low, 1);
% figure;
% stem(impulseResponse);
% title('FIR高通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 对音频信号进行高通滤波处理
% filteredAudioData_high = filter(b_high, 1, audioData);
% % 绘制高通滤波后的音频波形图
% figure;
% plot(t, filteredAudioData_high);
% xlabel('时间 (秒)');
% ylabel('幅度)');
% title('高通滤波后的音频波形图');
% % 对高通滤波后的音频信号进行FFT变换
% Y_high = fft(filteredAudioData_high);
% % 计算双边频谱 P2_high
% P2_high = abs(Y_high/length(filteredAudioData_high));
% % 计算单边频谱 P1_high
% P1_high = P2_high(1:length(P2_high)/2+1);
% P1_high(2:end-1) = 2*P1_high(2:end-1);
% % 绘制高通滤波后的音频频谱图
% figure;
% plot(f, P1_high);
% xlabel('频率 (Hz)');
% ylabel('幅度)');
% title('FIR高通滤波后的音频频谱图');
% %% butter高通
% % 计算滤波器阶数
% [N, Wn] = buttord(Wp, Ws, Rp, Rs);
% N = min(N, 4); % 限制滤波器阶数
% % 使用SOS表示设计滤波器
% [sos, g] = butter(N, Wn, 'high');
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('巴特沃斯高通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('巴特沃斯高通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('巴特沃斯高通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('巴特沃斯高通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_high = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% t = (0:length(audioData)-1)/fs; % 时间向量
% figure;
% plot(t, filteredAudioData_high);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Butterworth高通滤波后的音频波形图');
% % 频域分析
% Y_high = fft(filteredAudioData_high);
% P2_high = abs(Y_high/length(filteredAudioData_high));
% P1_high = P2_high(1:length(P2_high)/2+1);
% P1_high(2:end-1) = 2*P1_high(2:end-1);
% % 绘制频谱图
% figure;
% plot(f, P1_high);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Butterworth高通滤波后的音频频谱图');
% %% cheby1高通
% % 计算滤波器阶数N
% [N, Wn] = cheb1ord(Wp, Ws, Rp, Rs);
% N = min(N, 8); % 限制滤波器阶数
% % 设计Chebyshev I低通滤波器
% [sos, g] = cheby1(N, Rp, Wn, 'high');
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('cheby1高通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('cheby1高通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('cheby1高通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('cheby1高通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_high = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% figure;
% plot(t, filteredAudioData_high);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Chebyshev Type I高通滤波后的音频波形图');
% % 频域分析
% Y_high = fft(filteredAudioData_high);
% P2_high = abs(Y_high/length(filteredAudioData_high));
% P1_high = P2_high(1:length(P2_high)/2+1);
% P1_high(2:end-1) = 2*P1_high(2:end-1);
% % 绘制频谱图
% figure;
% plot(f, P1_high);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Chebyshev Type I高通滤波后的音频频谱图');
% %% ellip高通
% % 使用ellipord计算最小阶数
% [N, Wn] = ellipord(Wp, Ws, Rp, Rs);
% N = min(N, 4); % 限制滤波器阶数
% % 使用ellip函数设计滤波器
% [sos, g] = ellip(N, Rp, Rs, Wn, 'high'); % 使用SOS表示设计滤波器
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('ellip高通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('ellip高通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('ellip高通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('ellip高通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_high = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% figure;
% plot(t, filteredAudioData_high);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Elliptic高通滤波后的音频波形图');
% % 频域分析
% Y_high = fft(filteredAudioData_high);
% P2_high = abs(Y_low / length(filteredAudioData_high));
% P1_high = P2_high(1:length(P2_high) / 2 + 1);
% P1_high(2:end-1) = 2 * P1_high(2:end-1);
% % 绘制频谱图
% f = (0:length(P1_high)-1) * sampleRate / length(P1_high); % 频率向量
% figure;
% plot(f, P1_high);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Elliptic高通滤波后的音频频谱图');
% 
% %% FIR带通
% fp1 = 1200; % 下限通带截止频率（Hz）
% fp2 = 3000; % 上限通带截止频率（Hz）
% fs1 = 1000; % 下限阻带截止频率（Hz）
% fs2 = 3200; % 上限阻带截止频率（Hz）
% % 归一化频率
% Wp = [2*fp1/sampleRate 2*fp2/sampleRate];
% Ws = [2*fs1/sampleRate 2*fs2/sampleRate];
% % 计算过渡带宽
% DeltaW1 = Ws(2) - Wp(1);
% DeltaW2 = Wp(2) - Ws(1);
% DeltaW = min(DeltaW1, DeltaW2); % 使用min确保DeltaW为正数
% % 使用Hamming窗的经验公式计算滤波器阶数N
% N_band = ceil((6.6*pi) / DeltaW);
% % 设计FIR带通滤波器
% b_band = fir1(N_band, Wp, 'bandpass', hamming(N_band+1));
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(b_low, 1, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('FIR带通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('FIR带通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(b_low, 1, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('FIR带通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(b_low, 1);
% figure;
% stem(impulseResponse);
% title('FIR带通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 对音频信号进行带通滤波处理
% filteredAudioData_band = filter(b_band, 1, audioData);
% % 绘制带通滤波后的音频波形图
% figure;
% plot(t, filteredAudioData_band);
% xlabel('时间 (秒)');
% ylabel('幅度)');
% title('FIR带通滤波后的音频波形图');
% % 对带通滤波后的音频信号进行FFT变换
% Y_band = fft(filteredAudioData_band);
% % 计算双边频谱 P2_band
% P2_band = abs(Y_band/length(filteredAudioData_band));
% % 计算单边频谱 P1_band
% P1_band = P2_band(1:length(P2_band)/2+1);
% P1_band(2:end-1) = 2*P1_band(2:end-1);
% % 绘制带通滤波后的音频频谱图
% figure;
% plot(f, P1_band);
% xlabel('频率 (Hz)');
% ylabel('幅度)');
% title('FIR带通滤波后的音频频谱图');
% %% butter带通
% % 计算滤波器阶数
% [N, Wn] = buttord(Wp, Ws, Rp, Rs);
% N = min(N, 8); % 限制滤波器阶数
% % 使用SOS表示设计滤波器
% [sos, g] = butter(N, Wn, 'bandpass');
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('巴特沃斯带通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('巴特沃斯带通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('巴特沃斯带通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('巴特沃斯带通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_band = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% t = (0:length(audioData)-1)/fs; % 时间向量
% figure;
% plot(t, filteredAudioData_band);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Butterworth带通滤波后的音频波形图');
% % 频域分析
% Y_band = fft(filteredAudioData_band);
% P2_band = abs(Y_band/length(filteredAudioData_band));
% P1_band = P2_low(1:length(P2_band)/2+1);
% P1_band(2:end-1) = 2*P1_band(2:end-1);
% % 绘制频谱图
% figure;
% plot(f, P1_band);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Butterworth带通滤波后的音频频谱图');
% %% cheby1带通
% % 计算滤波器阶数N和截止频率Wn
% [N, Wn] = cheb1ord(Wp, Ws, Rp, Rs);
% N = min(N, 4); % 限制滤波器阶数
% % 设计Chebyshev Type I带通滤波器
% [sos, g] = cheby1(N, Rp, Wn, 'bandpass');
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('cheby1带通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('cheby1带通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('cheby1带通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('cheby1带通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_band = filtfilt(sos, g, audioData);
% figure;
% plot(t, filteredAudioData_band);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Chebyshev Type I带通滤波后的音频波形图');
% % 频域分析
% Y_band = fft(filteredAudioData_band);
% P2_band = abs(Y_band/length(filteredAudioData_band));
% P1_band = P2_band(1:length(P2_band)/2+1);
% P1_band(2:end-1) = 2*P1_band(2:end-1);
% figure;
% plot(f, P1_band);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Chebyshev Type I带通滤波后的音频频谱图');
% %% ellip带通
% % 使用ellipord计算最小阶数
% [N, Wn] = ellipord(Wp, Ws, Rp, Rs);
% N = min(N, 4); % 限制滤波器阶数
% % 使用ellip函数设计带通滤波器
% [sos, g] = ellip(N, Rp, Rs, Wn, 'bandpass'); % 使用SOS表示设计滤波器
% 
% %% 评估
% % 计算频率响应
% [H, f0] = freqz(sos, g, 1024, sampleRate);
% % 绘制幅度响应
% figure;
% subplot(2, 1, 1);
% plot(f0, abs(H));
% title('ellip带通滤波器的幅度响应');
% xlabel('频率 (Hz)');
% ylabel('幅度');
% % 绘制相位响应
% subplot(2, 1, 2);
% plot(f0, angle(H));
% title('ellip带通滤波器的相位响应');
% xlabel('频率 (Hz)');
% ylabel('相位 (弧度)');
% % 绘制群时延
% groupDelay = grpdelay(sos, g, 1024, sampleRate);
% figure;
% plot(f0, groupDelay);
% title('ellip带通滤波器的群时延');
% xlabel('频率 (Hz)');
% ylabel('群时延 (样本)');
% % 绘制滤波器的冲激响应
% impulseResponse = impz(sos, g);
% figure;
% stem(impulseResponse);
% title('ellip带通滤波器的冲激响应');
% xlabel('样本');
% ylabel('幅度');
% 
% % 使用filtfilt进行滤波
% filteredAudioData_bandpass = filtfilt(sos, g, audioData);
% % 绘制滤波后的波形
% t = (0:length(audioData)-1) / sampleRate; % 时间向量
% figure;
% plot(t, filteredAudioData_bandpass);
% xlabel('时间 (秒)');
% ylabel('幅度');
% title('Elliptic带通滤波后的音频波形图');
% % 频域分析
% Y_bandpass = fft(filteredAudioData_bandpass);
% P2_bandpass = abs(Y_bandpass / length(filteredAudioData_bandpass));
% P1_bandpass = P2_bandpass(1:length(P2_bandpass) / 2 + 1);
% P1_bandpass(2:end-1) = 2 * P1_bandpass(2:end-1);
% % 绘制频谱图
% f = (0:length(P1_bandpass)-1) * sampleRate / length(P1_bandpass); % 频率向量
% figure;
% plot(f, P1_bandpass);
% xlabel('频率 (Hz)');
% ylabel('幅度');
% title('Elliptic带通滤波后的音频频谱图');
